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Abstract. We consider a ID periodic atomistic model, for which we formulate and analyze 
an adaptive variant of a quasicontinuum method. We establish a posteriori error estimates for 
the energy norm and for the energy, based on a posteriori residual and stability estimates. 
We formulate adaptive mesh refinement algorithms based on these error estimators. Our 
numerical experiments indicate optimal convergence rates of these algorithms. 



1. Introduction 

Quasicontinuum (QC) methods, or more generally, atomistic-to-continuum coupling (a/c) 
methods, are a class of multiscale methods for coupling an atomistic model of a solid with a 
continuum model. These methods have been widely employed in atomistic simulations, where 
a fully atomistic model would result in prohibitive computational cost but is required in 
certain regions of interest (typically neighbourhoods of crystal defects) to achieve the desired 
degree of accuracy. A continuum model is applied to reduce the cost of the computation of 
the elastic far-fields. We refer to the recent review articles [L3], El] for introductions to QC 
methods. 

In this paper we present an a posteriori error analysis of a simplified variant of the energy- 
based QC approximation proposed in [25] . 

Considerable effort has been devoted to the a priori error analysis of QC methods [TTT [T2l 
[221 El El ESI [2D 123 ESI ESI ED] • A posteriori error control has received comparatively little 
attention: Arndt and Luskin use a goal-oriented approach for the a posteriori error estimates 
for the QC approximation of a Frenkel-Kontorova model [21 [31 E] . The estimates are used 
to optimize the choice of the atomistic region and the finite element mesh in the continuum 
region. Prudhomme et al. [21] also use the goal-oriented approach to provide a posteriori 
error control for the original energy-based QC approximation [17], which is inconsistent. Both 
of these approaches require the use of the solutions of dual problems. In [22J a posteriori 
error bounds for an energy norm are derived, using a similar approach as in the present work. 
However, the QC method analyzed in [22J does not contain an approximation of the stored 
energy and its 2D/3D version is therefore not a computationally efficient method. 

In the present work, we use the residual-based approach [281 Chapter 2], which is well 
established in the finite element approximation of partial differential equations, to derive a 
posteriori error bounds for the energy-norm and for the energy itself. What distinguishes 
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our setting from the classical one, are two particular features: 1. the model approximation 
("variational crime") is fundamentally different than quadrature approximations; 2. it is 
possible to adapt the mesh to a degree that the residual vanishes exactly in certain regions of 
interest (the atomistic region). Our work extends [IB], which considers a simplified setting. 

We do not tackle the question of model and mesh adaptivity for defect nucleation, but 
focus only on automatically choosing the size of the atomistic region and the finite element 
mesh in the continuum region. In many applications, a defect is inserted into the crystal 
before the computation, or it is known a priori where a defect will nucleate. An extension 
of our work to include defect nucleation would be valuable, but cannot be pursued in the 
simplified ID setting we are considering here. 

1.1. Outline. In Section[2]we introduce a ID atomistic model problem, which mimics the be- 
haviour of crystal defects in 2D/3D. Moreover, we review the construction of a QC method to 
efficiently approximate its solutions, and introduce the notation that will be used throughout 
the paper. 

In Section [3j we derive a residual estimate for the QC method in a discrete negative 
Sobolev norm. In Section [4j we present an a posteriori stability result. In Section [5j we 
combine the residual estimate and the stability result to give a posteriori error estimates for 
the deformation gradient and for the energy. 

In Section |6j we describe three mesh refinement algorithms based on our a posteriori 
error analysis and on previous a priori error estimates, and present a numerical example to 
illustrate the performance of these algorithms. 



2. Model Problem and QC Approximation 

2.1. Atomistic Model. Following previous works [7J [181 123] we formulate a model problem 
in a discrete periodic domain containing 2N atoms, where N G N. Let F > denote a 
macroscopic stretch and e = 1/(2N) the lattice spacing, both of which we fix throughout. 
We define the displacement and deformation spaces, respectively, by 

U £ : = {« G 1R Z : U£ +2 n = Ue, and u = 0}, and 
y £ := {y G R z : y t+2N = F + y t , and y = 0}. 

In particular, we observe that y G y e if and only if y = Fx + u for some u 6W £ . 
The stored energy (per period) of an admissible deformation y G y £ is given by 

N N 

S a {y):=e ^ <f){y' e ) + e M-i+Vi)' 
e=-N+i e=-N+i 

where y' e := e~ 1 (^ — ye-i) and we note that y\_ x + y\ = s^ 1 (ye — ye-2) is a second-neighbour 
difference, and where G C 3 (0,+oo) is a Lennard- Jones type interaction potential. We 
assume throughout that there exists r* > such that is convex in (0, r*) and concave in 
(r#, +00). 
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Given a periodic dead load / G C (M), f(x + 1) = f(x), we define the external energy (per 
period) by 

N 

-(f,u) e :=-e f« u *> 
e=-N+i 

where u = y — Fx . Thus, the total energy (per period) under a deformation y G y e is given 
by 

E*(y):=£ A (y)-(f,y-Fx) e . 

We wish to compute 

y a G argminE a (y e ), (2.2) 
where argmin denotes the set of local minimizers. 

Remark 1. 1. Since the internal energy is translation invariant, our choice uq = (instead 
of the more common constraint Ylf=-N+i u £ = 0) does not alter the problem but simplifies 
the treatment of external forces in Sections 13.21 and 13. 41 



2. The external energy (2.1) can be thought of as the linearisation of a nonlinear external 
energy about the state Fx. 

3. We have included the external forces primarily to render the problem non-trivial. In 
realistic applications in 2D/3D, defects or forces applied on a boundary are the cause of 
deformation of the crystal, however, in ID such complex behaviour cannot be described 
directly. □ 

2.2. Notation for lattice functions. Throughout, we identify a lattice function v G M. z 
with its continuous and piecewise affine interpolant with nodal values v(e£) = V£. Vice- versa, 
if g G C°(R) then we always identify it with its vectors of nodal values (gi)eez — (g( £ ^))tez- 
Let D be a subset of Z. For a vector v G M z , we define 



maxe & T)\ve\, p = oo. 



\ v \\g<p) :z 



If the label T> is omitted, then we understand this to mean T> = {—N + 1, . . . , N}. 

We define the first discrete derivatives v' £ := (v£ — ve-i)/s. It is then straightforward to see 
that, if at the same time we identify v' with the pointwise derivative, then \\v '||lp(-i/2,i/2) = 
IMI«- 

We equip the space U e with the discrete Sobolev norm (recall that u = 0) 

IMIw 1 - 2 := \W I U 2 (-1/2,1/2) for v e u 6 - 

The norm on the dual (U £ )* is defined by 

||T|| w -i, 2 := sup T[v}. 

IMI M 1,2=1 
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2.3. Finite element notation. To construct the QC approximation in the next section, we 
first define some convenient notation. Let Q := [—1/2,1/2] denote the computational cell. 
We choose an atomistic region f2 a C Q, where atomistic accuracy is required, and we define 
the continuum region by Q c := 0\O a . We will also use the periodic extension of Q c , denoted 
by n* := (n c + Z). 

We assume throughout that fl a is an open interval (L a , R a ) with —1/2 < L a < i? a < 1/2. 



All our results (with the exception of Section 3.4) can be extended without difficulty to the 



case when f2 a consists of a finite union of open intervals. 

Let {T^} fceZ be a partition of R into closed intervals with Tj? = [Xk-n x k)i where x\ > x\_ x 
are the nodes of the partition. We assume, without loss of generality, that x\ is the left-most 
node and x h K the right-most nodes in the interval (—1/2, 1/2]. 

The length of an element is denoted by h k := \Tjj\ := x\ — x\_ x . The space of continuous 
piecewise affine functions with respect to the partition T h is denoted by V\(T h )- 

We assume throughout that the partition T h has the following properties: 

(Tl) T h is periodic: there exists K G N such that x k+K = 1 + x\ for all k G Z. 
(T2) T h has atomistic resolution in fl a : Q a n eZ = Q a fl {x k }kez- 

(T3) The a/c interface points are finite element nodes: <9f2 a C {xfykez- In particular, each 

element belongs entirely to either the atomistic or continuum region. 
(T4) If T£ C Q# then \T£\ = h k > 2e. 

Property (T4) is not strictly required, but simplifies the analysis and is not a significant 
restriction. Note also that we have not required (except in the atomistic region) that finite 
element nodes must be positioned on atomic sites. Although not necessary in ID, it somewhat 
simplifies mesh generation, and to some extend mimics the fact that element edges or faces 
in 2D/3D cannot normally be aligned with the underlying crystal lattice. 

The finite element displacement and deformation spaces are defined, respectively, by 

U h := {u h G Vi{T h ) : u h {x + 1) = u h {x) and u h {0) = 0}, and (2.3) 

y h ■= {y h G 7>i(7*) -Vh-Fxe U h }. (2.4) 

For g G C°(R), we define the interpolation operator Ih ■ C°(R) — >■ Vi(T h ) by 

(I h g)(x h k ) := g(x h k ) \/k G Z. (2.5) 

We note that I h : W ->• U h . 

For future reference, we also define the micro-elements Tf := ((£ — l)e,£e) for i G Z. 
Analogously, we define I e to be the nodal interpolant with respect to the atomistic grid. We 
will require this interpolant since the mesh nodes {x k } k ez do not necessarily coincide with 
lattice sites. 



2.4. QC Approximation. The QC approximation we analyze in this paper is the ID variant 
of the ACC method described in [26]. (An earlier variant of the idea was described in [18] 
and a similar construction in [10] . We focus on the formulation proposed in [26J since it can 
be readily generalised to 2D.) 

The idea of the ACC method is based on the splitting of interaction bonds. A bond is an 
open interval b = (£e,(£ + j)e) for f 6 Z and j G {1,2} (since we consider only first and 
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second neighbour interactions). Since our computational domain is (0, 1] the set of bonds 
over which the atomistic energy is defined is given by 

B := { {£e, {£ + j)e) : j = 1, 2; 1= -N + 1, . . . , N}. (2.6) 

For each bond b = (£e, (£ + j)e) we define r b := j. 

For any open interval u = (L u , R u ) (e.g., for a bond) with length \u>\ := R w — L w > we 
define the finite difference operator 

DuV : = vJM^M for v e C °(R). (2.7) 

Note that, with this notation, we can rewrite £a,(y) = £ Ylib&B 4>{ r bDby) ■ 

For each bond b and deformation field y G W /1,00 (1R) we define its atomistic and continuum 
energy contributions to the stored energy, respectively, by 

a>b(v) bC !^ a ^{nDbnn^y) and c b (y) := j- [ ^ <j)(r b y'(x)) dx. 



# 



bna 

If \b fl O a | = then D^^y is ill-defined; in that case we define a b = 0. The QC energy (the 
ACC energy of [2S]) of a deformation field y G W 1 ' 00 (R) is now defined by 

£ qc (y) :=eY,[ab(y)+c b (y)]. (2.8) 

beB 

The key property of this definition is that it satisfies the patch test [2S1 Section 3.3], 

£' & (Fx)[v £ } = £' qc (Fx)[v h ] = for all v £ G IA £ and Vh G Li h , VF > 0. (2.9) 
The external energy (per period) is given by 

-(f,u h ) h :=- [ I h (fu h )dx. (2.10) 



Note that, in the atomistic region, this reduces to the same form as (2.1). The total energy 
(per period) of a deformation y^ G y h is then given by 

E qc (y h ) ■= £ qc (yh) - (/, yh - Fx) h . 

In the QC approximation we seek 

y qc e&rgmm E qc (y h ). (2.11) 



Remark 2. It is initially not obvious why the formulation (2.8 ) should reduce the complexity 
of the computation of y qc over that of y a , since £ qc is still written as a summation over all 
bonds. However, one can readily check (see [2S] for the details), that 



£qc{yh) = a ^ yh "> + / W (y'h) dx ' 
beB ^ Qc 



where W(r) := 0(r) + 0(2r) is the Cauchy-Born stored energy function. This formulation 
requires only a sum over all bonds within the atomistic region, and a standard finite element 
assembly procedure to compute the energy contribution from the continuum region. Thus, 
the evaluation of £ qc has only a computational complexity proportional to #7" . 
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Analogous results hold in 2D [25]; a more involved assembly procedure is required to make 
the ACC method efficient in 3D EH] . □ 



3. Residual Analysis 



3.1. The residual. Let y & be a solution of the atomistic problem (2.2). It is straightforward 



to see that, if y' a > 0, then £ a is differentiable at y a and hence the first order optimality 



condition for (2.2) is satisfied: 

£'M)[v] = (f,v) £ \fveu\ (3.1) 

where [v] = e ^ 0'(r 6 D 6 y a )ryDj,T;. (3.2) 



beB 



Similarly, if y qc is a solution of the QC problem (2.11 ) with y' > on [—1/2, 1/2], then it 
satisfies the corresponding first order optimality condition 

£' qc (y qc )iv h ] = (f,v h ) h Wv h eu\ (3.3) 



where £' qc {y h )[vh\ = Y (a b (y h )[v h ] + c' b (y h )[v h }) (3.4) 

beB 

= Y\ bn ft a \4>'(r b D br)Ua y h )D bnn& v h + Y # ^ \ r dh( x )) v 'h dx - 



beB beB 



Let y qc G y h be a solution of (2.11), then we define its residual R G (U £ )* by 

fl[«] := E' a (y qc )[v}, (3.5) 



Using (3.3) we can rewrite this as 

R{v] = E' a (y qc )[v} - E' qc (y qc )[I h v] 

= {S'M[v} - £' qc (y qc )[I h v}} - {{f,v) £ -(f,I h v) h } 

= : -RjntM + -RextH, 

where we call i?j nt the internal residual and R ext the external residual. We will bound them 
separately in the next two sections. 

3.2. Estimate for the internal residual. In this section, we analyze the internal residual 
R m t- To state the main theorem, we define the index set of all nodes in the continuum region 
and a/c interface (recall that Q c is closed), 

/C c := {k G {1,...,K} : x h k G fi c }. (3.6) 

Theorem 1. Let y h G y h such that y' h > and Ri nt [v] = £' & (y h )[v\ - £' qc (y h )[I h v], then 

||-Rmt|| M -i,a < ( 3 Y Vl) 2 ='■ ViVh), where (3.7) 
ke/Cc 

rfk '■= Y (W^bDbVh) ~ 4>'(r b D bnn!L y h )\\ 2 L2(bnna) + \\(f>'(r b D b y h ) - <p'{r b y' h ) ||* 2(6nn # } ) • 



beB 



x%eint(b) 
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Remark 3. 1. The expressions for rjk are reminiscent of the flux (or stress) jump terms that 
occur in the classical residual error analysis for partial differential equations. The origin in 
our case, is somewhat different however, and results only from the model approximation and 
not the finite element coarsening. 

2. With some additional work, the form of the estimates pk can be turned into element 
contributions and further simplified by computing more explicit representations. □ 



Proof. Let Vh '■= Ih v - From (3.2) and (3.4) we obtain 

Rmt[v] = y2\ \b\4>'(r b D b y h )D b v - \b H ft a |0' '(r b D bnUa y h )D bnUa v h - / (f)' (r b y' h )v' h dx\. 

(3.8) 

We subtract and add the terms 

y~] 1 6 n ^l a \(/)'(r b D bn u a yh)Dbnn a v and ^ / <P'( r b!/h( x )) v ' dx 
beB beB J bnn c 

to split -Ri n t into three components, 

= ^2 \ \b\<fi'{nD b y h )D b v -\b<l Q a \<p'{r b D bnnii y h )D bnnii v - / <f>' (r b y' h (x))v ' dx 
beB I Jbnn* 



+ ^2\bn ^WinD^^yh) (D bnQa v - D bnn& v h ) 
beB 

+ V / J)'(r b y' h (x))(v' -v' h ) dx 
Jbnn* 



ibnnjf 

=: RM + R 2 [v] + R 3 [v\. (3.9) 

We will show that R 2 = -R3 = and estimate 

For i? 2 this follows simply from the fact that v h = v in f2 a and hence D bn n a v = D bn Q & Vh 
for all 6 G B such that |6 n Q a \ > 0. 

For i? 3 , following the computations in [26], Section 3.2], we obtain 



R 3 [ V ]= [ W'(y' h )(v'-v' h )dx, 
Jn r 



where W{r) = 0(r) + 0(2r). Since v(x%) = Vh(x^) for all k G Z and since W(y^) is constant 
on each element it follows that -R3 = 0. 

Finally, we turn to the analysis of R\ = Y^beB R\i wnere we define 

R\[v\ := \b\(j)' (r b D b y h )D b v - |& n Q a |<//(r 6 Ann a ^)A>nn a ^ - / (j)' {r b y' h {x))v' dx. 

Jbnn* 
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Using the fact that \oj\DuV = f v' dx we obtain 

R b 1 [v]= (f)'(r b D b y h )v' dx - <p' (r b D br] ^y h )v' dx - / (jj ' (r b y' h (x))v' dx 
Jb Ann a Jbnn* 

[(j)'(r b D b y h ) - (f)'(r b D bnna y h )] v dx + [0 \r b D b y h ) - (j)' (r b y' h )]v' dx. 
'bnn a Jbnnf 

If b C Vl a , then b n fi a = b and \b n f2 c | = 0, and hence .Rj = 0. Similarly, if 6 C C fif , 
then Z?b2//i = 2//Jt' 1 an d 1^ H O a | = and hence i?J = 0. Thus, we observe that only bonds 
crossing continuum element boundaries, or the atomistic/continuum interface, contribute to 
the residual. These are precisely the points contained in /C c . In particular, we obtain 

RM = E E ^H> (3.io) 

fceKc beB 

x^eint(fe) 

where we used the fact that no bond can cross more than one point x\ G /C c due to our 
assumption that all elements have at least length 2e. 

From the definition of R\, and applying two Cauchy-Schwarz inequalities, it is straight- 
forward to estimate 

|2 

lz, 2 (ftnn a ) 



Ri[v}\ < ( \\4>'(r b D b y h ) - (j)' {r b D bm& y h )\ 

+ \\<t>' {r b D b y h ) - (j) (r&J/A)|| £ a (6nn # ) J IMU 2 (&), 



1/2 

J|L 2 (&nnf ) J 

and after applying another Cauchy-Schwarz inequality 



E ^m|<%( E 

beB beB 

x£eint(6) x£eint(&) 



I 'II 2 

\ V llL 2 (6) 



1/2 



where % is defined in (3.7). 

Combing our foregoing estimates, we arrive at 

1/2 



Rim < e>( E 



i 

l v llL 2 (b) 



x£eint(6) 



<'E^f(E E 



x%eint(b) 



and we are only left to estimate the sums involving the test function. 

To that end, we simply note that, due to (T4), for any fixed point x G (—1/2,1/2], the 
maximal number of bonds appearing in the sum on the left-hand side below and crossing x 
is three; hence, 

E E W v '\\h(b)<3\W\\h- 
keK. c beB 

x£eint(6) 

This concludes the proof. □ 
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3.3. Estimate of the external residual. We now turn to the estimate of the residual of 
the external energy, which was defined as 

R ext [v] = (f,v) £ - (f,I h v) h = [ [I e (fv) - h(fhv)] dx. (3.11) 

Jn 

We outline the key points of the argument for estimating -R ex t, before stating the result. 
We define a slightly extended continuum region, 

&c--=\J{T!: |T/ n Q c | >0}; 
then h(fv) = h{fhv) in f2 a \ Q c , and therefore 
R cxt [v}= f [l E (fv)-I h (fhv)]dx 

= L IWv) - fv] dx + I [fv - fl h v] dx + ! [fl h v - h{fhv)] dx 
=:R 1 [v] + R 2 [v]+R 3 [v]. (3.12) 

The three terms can be estimated using standard interpolation error results, hence we only 
give a brief outline of the proof of the resulting bound. 

Proposition 2. Let f G C 2 {Vt c ), then 

<V f + V q , (3-13) 

where the error due to external forces ly and the "quadrature error" rj q are, respectively, 
defined as follows: for * G {/, q}, 



k&{l,...,K} 
J\2 ._ h 2 k I, ,„2 



where {r,^ 2 == ^ll/ll^e*)) 



and (4)' := (e 4 + hi)\\f\\* h + ^llf " 2 



Remark 4- 1- Note that there is an error contribution from the atomistic region, due to the 
fact that in the elements touching the a/c interface, the "quadrature" approximation of the 
external forces is not exact. For the purpose of mesh refinement, we count this error towards 
the neighbouring elements in the continuum region. 

2. An alternative residual estimate that does not use / G C 2 (Cl c ), but only the discrete 
setting, is presented in [29]. This requires a much more involved argument. □ 



Proof. From (3.12) we obtain 



R ex t[v]=Ri[v]+R 2 [v]+R 3 [v]. 
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Applying standard interpolation error estimates (see, e.g., |5J 0]) on elements T|,T^, we 
obtain 

/ [h(fv) - fv] dx < l\\e 2 f"\\ L 2 iT!) \\v\\ L 2 {T!) + |||£ 2 /'IU 2 (r/) II^H^^i), 

l[f V ~ f J h v ] dx < $\\ h kf\\tf(Tj:)\\ v '\\l?(T£)> and 

[h{fv) - fl h v] dx < l\\h 2 J" W^^WhvW^^ + iWhlf'W^^Wlhv'W^^hy 



Summing over all elements, applying the Cauchy-Schwarz inequality, and defining h(x) := hk 
for x G Tj}, 

Ri[v] < iWf'Wvfa) IMU 2 + y 11/11^(0:) W v 'Wl 2 ' 

RM < iWhfW&fa) \WW&> and 

Rs[v] < l\\h 2 f"\\ L2{Clc) \\hv\\ L * + l\\h 2 f\\ L ^ c) \\I h v'\\v>. 

We now use the estimates (which exploit the fact that l^v' is the L 2 -orthogonal projection 
of v' onto piecewise constants) 

IMIl2 < JIKIU*, \\hv'\\u> < \w\\&, 

and \\I h v\\ L 2 < l\\I h v'\\ L 2 < ±\\v'\\ L 2, 

to deduce that 

< Bhf\\ LH nJv'\\L* + (^||(£ 2 + h 2 )f"\\ 2 LH n c) + 11^ + h 2 )ni^ /2 \W\W- 
The result follows by splitting the norms inside the brackets over elements. □ 

3.4. External residual estimate for singular forces. In our numerical experiments in 
Section [6] we shall employ an external force that behaves like ~ N _1 near x = 

(we use the "singularity" in the force to mimic a defect). Let us suppose that we also have 
| /'(#)! ~ \x\~ 2 and \f"(x)\ ~ \x\~ 3 near the origin. We now give a formal motivation why 
the quadrature estimates employed in Proposition [2] are inadequate in this situation. 
Applying the quadrature estimates to such a force field, we obtain 

4 - h 2 \\r\\ L , in) + h 2 \\r\\ L2{Tkh) » hT\x h k \-" 2 + hT\x h k \-^ 2 . 

We notice that, for T k near the origin, H/"!^^^ ~ Htl -1 !!/' Hz^t' 1 )- Moreover, the quadra- 
ture estimate is 0(1) and cannot be controlled. By contrast, 

r,f ~ A, 3/ 2 1 I — 1 
'Ik ~ a k l x fcl ' 

from which we conclude that hlWf'W^^h^ is dominated by ry{, but that r][ is itself dominated 

by h 2 \\f"\\ L 2 {T hy 
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The origin of this undesirable effect is the (ab-)use of the Poincare inequality in the proof 
of Proposition [2] In the remainder of this section, we shall remedy the situation by replacing 
the standard Poincare inequality with a weighted variant. This approach is inspired by [20] . 



Lemma 3. Let £l c be defined as in Section 3.3, and let w(x) := xlog 2 (x) ; then 



Proof. We begin by noting that 

\v(x)\ < {x^^Wv'Wtf for all iGfi. 

Hence, we can estimate 

,-1/2 rl/2 1 

/ {w^vix^ 2 dx < Wv'Wjpp 1/2) / 2 dr. 

Jr. Jr. x log (2x) 

Since (log _1 (x))' = — (x log 2 (x)) _1 we obtain 

yl/2 

Applying an analogous argument in the left half of the domain, we obtain the stated estimate. 

□ 

We now apply this estimate to obtain an alternative external residual estimate. 

Proposition 4. Let f 6 C 2 (Vt c ), then 

Pextllw-M <v f + V q , (3-14) 
where is defined in (3.13) and fj q is defined as follows: 



m 2 --= £ 



g\2 
k) > 



fee{i,...,-fsr} 

^ere (r^) 2 := (e 4 + ^ll/'ll^) + ^IKlkp*)' 
and where w is defined in Lemma [3| 



Proof. We again use the splitting (3.12) to obtain 

#extM = Ri[v]+R 2 [v] + R 3 [v). 

The residual term R^o\ is estimated in the same way as in the proof of Proposition [2j and 
gives rise to the term in the estimate. 

We show only the modified estimate for i^fi*], since the estimate for R\[v\ is analogous. 
Applying again a standard interpolation error estimate, we obtain 

\Rs[v]\ < l\\h 2 (fI h v)"\\ L ^ c) < l\\h 2 f'I h v'\\ Ll ^ c) + l\\h 2 f"I h v\\ Ll{flc y (3.15) 

The term \\\h 2 f Ihv'W^^^ can be treated in the same way as in the proof of Proposition 
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To estimate the second term on the right-hand side of (3.15) we insert the weighting 



function w defined in Lemma [3] and then apply the weighted Poincare inequality: 

l\\h 2 f"I h v\\ Ll{Clc) = \\\(h 2 wf"){w- l I h v)\\ L1(Slc) 

<\\\h 2 wf"\\ L *(n c ) lk rl 4^|| L 2 ( ^ c) 

<4^2\\ h2w f"\\L2(Cl c ) WW*. 

By continuing to argue as in the proof of Proposition [2] we obtain the stated estimate. □ 

We can now revisit the issue of relative magnitude of the various contributions to the 
residual estimate for the case where \f(x) \ ~ ~ \ x \~ 2 an d \ f"i x ) \ ~ \x\~ 3 . Note 

that the effect of the weighting function is that \w(x)f"(x)\ ~ | log 2 (x)| \x\~ 2 which is now 
comparable to U P to a 1°S factor. 

More precisely, suppose that Tj? is near the origin, then we now obtain 

^-^ /2 i4r 3/2 + ^ /2 i4r 3/2 iog 2 (4). 

In particular, we observe that in the new external residual estimate, the quadrature error is 
dominated by the main error term , which is the same as in the standard estimate given in 



Proposition 3.3 Thus, in our numerical algorithms presented in Section^we will be justified 



in neglecting the effect of the quadrature errors. 

4. Stability 

Stability of the exact (i.e., the atomistic) model is the second key ingredient for deriving 
an a posteriori error bound. Our aim is to prove coercivity (or, positivity) of the atomistic 
Hessian at the QC solution y qc : 

KW)M > c a (y q c)IKIl! 2 e u, (4.1) 

for some constant c a (y qc ) > 0, where the Hessian operator of the atomistic model is given by 

TV N 

K{y)[v,v] = e £ <P"(y'iM\ 2 + z E ny'i + ve + iU + v' e+1 \ 2 . 

e=-N+i e=-N+i 

Following [5J [TS] we note that the 'non-local' Hessian terms \v\ + f£ +1 | 2 can be rewritten in 
terms of the 'local' terms \v' £ \ 2 and |^ +1 | 2 and a strain-gradient correction, 

|^ + ^ +1 | 2 = 2|^| 2 + 2|^ +1 | 2 -e 2 |^| 2 . 

Using this formula, we can rewrite the Hessian in the form 

JV N 

K{y)[vM=e £ A e \v>\ 2 + e £ Bi \v>tf, 
e=-N+i e=-N+i 

where 

A e (y) := cf>"( y ' t ) + 2^(y' i _ 1 + y' e ) + 2^(y' e + y' i+1 ) (4.2) 
Be(y) ■= - ftXy'e + y'e+i)- 



Recall our assumption in £2.1 that <p is convex in (0,7"*) and concave in (r*,+oo). For 



typical interactions such as Lennard- Jones or Morse potentials one generally observes that 
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y'e > r */2, hence we shall assume this throughout. As a result of this assumption, and the 
properties of (f), we have 5 f > Vf 6 Z. 

As an immediate consequence we obtain the following lemma, which gives sufficient con- 
ditions for the stability of the atomistic hessian evaluated at the QC solution. 



Proposition 5. Let y qc G y e satisfy mim^y^)^ > r*/2 ; th 



en 



Vf G U, where A*(y qc ) :— min A f (y, 



„ qc/j 

1,...JV 



and the coefficients Ai(y qc ) are defined in (4.2). 



Remark 5. Since the minimum in the definition of A* is taken over 2N lattice sites, it 
appears at first glance that A* is expensive to evaluate. However, exploiting the fact that y qc 
is piecewise affine, one can evaluate A* in O(K) operations: 

Case i: If dist(e(£ — 1/2), {x^}) < |e then we evaluate Ai(y qc ) using (4.2). There are 
0(K) lattice sites of this type. 

Case n: If dist(e(£ - 1/2), {x h k }) > |e then 

^(2/qc) = WJr^+VWJ^), 
that is, we only need to evaluate this formula once for each element. □ 

4.1. Estimates for the hessian. Before we present our main theorems, we state two use- 
ful auxiliary results: a local bound and a local Lipschitz bound on E". The proofs are 
straightforward and are therefore omitted. 

Lemma 6. Let y,z G y e such that min^ y' e > \i and min^ z' e > \x for some constant fi > 0, 
then 

\£"(y)[v,w]\ KCnWv'WpJw'Wp^ and 
\{E'l(y) - £"(z)}[v,w]\ <C Up \\y' - z'\\e r \\v'\\ e 2\\w'\\ e 2 V«,tw G 14, 
where C H := M 2 (/i) + 4M 2 (2 j u) and C Lip ■= M 3 (/x) + 8M 3 (2/x) and Mj(t) : = max s > t |0 (i) (s)|. 

5. A Posteriori Error Estimates 

5.1. A posteriori error estimate for the solution. We will assume the existence of an 
atomistic solution y & in a neighbourhood of y qc (cf. (5.1)), and estimate the error y & — 
I £ y qc - It is in principle possible to rigorously prove the existence of such a solution y a in a 
neighbourhoodo f y qc , following for example [221 [18] . however, this would require substantial 
additional technicalities. 



Theorem 7. Let y qc be a solution of the QC problem (2.11) with miae(y qc )' e > r*/2 and 
A*(y qc ) > 0, where A* is defined in the statement of Lemma^E, Suppose, further, that y & is 
a solution of the atomistic model (2.2) such that, for some t > 0, 

lbl-2/qclU°°(n) < r. (5.1) 
If t is sufficiently small, then we have the error estimate 

\\y' a -iey' qc \\L2(ti) < Ajj£)(v(v<tp) + v f + v q ), (5-2) 



14 CHRISTOPH ORTNER AND HAO WANG 



where rj is defined in (3.7) and rj* and r] q are defined in (3.13). 

Proof. Let e := y a — I £ y qc - We require that r < ^miny' qc , then by the mean value theorem 
we know that there exists 9 G conv{y a , I e y qc }, such that 

E'M[e,e] =E'M[e] - E' a (I £ y qc )[e] 

= E> qc (y qc )[I h e}-E' & (I £ y qc )[e] 

= (£' qc (y qc )[I h e] - £' a {I £ y qc )[e\) - ((f,I h e) h -(f,e) E ). 

In this proof we write £a.(I £ y qc ) to emphasize that we are comparing y' & with I £ y' qc . 
Applying Theorem [T] and Proposition [2j the two groups are respectively bounded by 

Kc(2/qc)[4e] - £' a (I E y qc )[e}\ < r]{y qc ) \\e'\\ L 2, and 

|</,4e) ft -(/,e> e | <(v f + V q )\W\\ L ^ 

and hence we obtain 

£MM< Mv*) + n , + ri q )\W\\»- ( 5 - 3 ) 

Next, we compute a lower bound on £ a (9)[e,e]. Using the Lipschitz estimate given in 
Lemma [6| Proposition [5] together with our assumption that y' qc > r*/2, and the a priori 



bound (5.1), we have 



£ a (0)[e,e\ >^'(y qc )[e,e] -C L ip||y , a -y / qc |U«.||e'|| La 
> (A^y qc )-C Up r)\\e'\\ 2 L2 . 

Hence, we if require that r < A*(y qc )/ (2Cu P ) (since Cu P is bounded as r decreases this is 
satisfied for r sufficiently small), then we arrive at 

iMVvWWb < {ri e (y v )+ri , + Ti 9 )\W\\»- 

Dividing through by || e'|| j^a yields the stated estimate. □ 

5.2. A posteriori error estimate for the energy. An important quantity of interest is 
the total energy of the system being approximated. In this section, we derive an a posteriori 
estimate for the energy difference E a (y a ) — E qc (y qc ). To that end we decompose the energy 
difference into 



\E a (y a ) - E qc (y qc ) \ = \E a (y a ) - E a (I £ y qc )\ + \£ a (I £ y qc ) - £ qc {y qc ) 

+ | (/, IsVqc ~ Fx) £ - (/, y qc - Fx) h \ 



(5.4) 



and analyze each component separately. 



For the first group on the right-hand side of (5.4), the result is standard: 



Lemma 8. Let y,z G y e such that min^ y' e > \x and min^ z\ > [i for some constant \x > 0, 
and y G argmin£' a (3^ e ), then 

\E a (y) - E a (z)\ < lC u \\y' - z'Wi*, (5.5) 
where Ch = Ch(aO was defined in Lemma^ 
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Proof. There exists 9 G conv{?/, z} such that 

EJiz) = E a (y) + E'M[z -y} + \E'i{e)[z -y,z-y]. 
Since E' a (y)[z — y] = and E" = E" applying Lemma [6] yields the stated result. 



□ 

For the second group on the right-hand side of ( ]5.4[ ), the estimate is obtained from a 
straightforward computation, using only the fact that the energy of a bond lying entirely 
inside an element is exact in the QC energy. The proof is omitted. Although the form of 
the error contributions fik looks complex at first glance, they are in fact straightforward to 
compute. 



Lemma 9. For y^ G y h and y' h > 0, we have 
\£*{yh) - £qc(yh) \ < KVh) ■= 22 /ifc ' where 



(5.6) 



£ 

x^eint(6) 



|bnQ a | 



1 

n Jbna* 



[(f)(r b D b y h ) - <p{r b y' h )] dx 



Finally, the third group on the right-hand side of (5.4) can be estimated similarly as the 
external residual, however since the "test function" Uh = yh — Fx is now known explicitly, 
some additional structure can be exploited. Note that the error due to quadrature is again 
of higher order. 



Lemma 10. Let y^ G y h and Uh '■= yh — Ex, then 

\(f,Uh)h - (f,I £ Uh)e\ < ti f (yh) + fi q (y h ), 



(5.7) 



where ^ (yh) is the error due to external forces, and n q (yh) the error due to quadrature. They 
are, respectively, defined as follows: 



l* f (Vh) 



where d 



f ._ 



& Where 



L^ k )\[u h \ x h\, and 

\\(f!eUh) H^i(Th) + ~t || (f u h) ||li(t' 1 )' 



where the second derivative (fI £ Uh)" is to be understood in the piecewise sense, = Tf if 
(£—l)e < x\ < is anduk = if no such (£— 1) G Z exists, and [u' h ] x h := u' h (x l l+)— u' h (x\— ). 



Remark 6. 1. It is straightforward to evaluate (possibly upper bounds of) the error 
contributions rj( and rfj, with 0(K) computational complexity. This is due to the fact that 
Uh is piecewise linear on Tj}, and I £ Uh = Uh except in the neighourhoods Uk of the element 
interfaces. 

2. For the purpose of mesh refinement, we will group the residual contribution of the 
elements touching the a/c interface, so that no error contributions is associated with the 
atomistic region, which cannot be further refined. □ 
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Proof. Similarly as in the external residual estimate we write the external energy difference 



as 



(f,IeUh)e - (f,U h ) h 



[le(fleUh) ~ Ih(fUh)] dx. 



Using I £ (fv) = Ih(flh,v) in Q \ Q c , we decompose this difference into 



[ie(fleUh) - h(fuh)] dx 



[ie(fleUh) -Ih(fUh)] dx 

[l £ (fl £ u h ) - fl £ u h ] dx+ [fhuh ~ fu h ] dx (5.8) 



+ / [fuh ~ h{fu h )] dx 
Jn c 



Unlike in the proof proof of Proposition [2j where Vh was unknown, we can use our explicit 
knowledge of Uh, to estimate the first and third groups on the right-hand side of (5.8) as 
follows: 



[l £ {fl £ u h ) - fl £ u h ] dx 
[l h (fu h ) - fu h ] dx 



<||£(/W|L 



L 1 ' 



and 



<IIt(WL 



(5.9) 



L 1 ' 



where the second derivatives are understood in a piecewise sense. 

To estimate the second group on the right-hand side of (5.8 ), we note that I £ Uh = Uh except 
near element boundaries. Upon defining Uk and [u' h ] Xk as in the statement of the result, we 
have 



/ 

J fir 



[fl £ u h - fu h ] dx = / [ftsUh - fu h ] dx 

k=l,...,K 

£ E 



< 



E S 

k=l,...,K 



L^(u k )\\U h hu' h \\ L(x ^ k) 



Li(oj k )\[U h \x k 



(5.10) 
(5.11) 



It is now straightforward to rearrange the various error contributions to obtain the stated 
result. □ 

Combining all the foregoing estimates yields an a posteriori error estimate for the energy. 
Theorem 11. Suppose that the conditions of Theorem^ are satisfied, then 



|J5 a (y a ) - E qc (y qc )\ < 



4C H 

^*(?/qc) : 



{V(y q c) 2 + (V f ) 2 + (V") 2 } + MZ/qc) + //(j/qc) + ^(j/qc), (5.12) 



where r](y qc ) is defined in fl3~7| ) ; r] f ,r] q in fl3.13[ ) ; fi(y qc ) in ( |5~6| ) and fi f (y qc ), fi q (y qc ) in ( |5~7| ). 
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Proof. The second term on the right-hand side of ( 5.4[ ) is estimated by ^(yn), which gives 



rise to the second term on the right-hand side of (5.12) (cf. Lemma Em. The third term on 



the right-hand side of (5.4) is estimated by fi*(yh) and fi Q (yh), which gives rise to the third 



and fourth terms on the right-hand side of (5.4) (cf. Lemma 10). 



Finally, using Lemma |8j the first term on the right-hand side of (5.4) can be bounded 
above by 

\E a (y & ) - E a (I £ y qc )\ < \C n \\y' a - hy'Jtf- 
Employing Theorem [7] we obtain 

\EM - E a (I £ y qc )\ < j^iviVvT + (v f f + (V q ) 2 )- □ 

6. Numerical Experiments 

In this section, we present numerical experiments to illustrate the results of our analysis, 
and highlight further features of our a posteriori error estimates. In particular, we will 
propose an adaptive mesh refinement algorithm, and show numerically that it achieves an 
optimal "convergence rate" in terms of the number of degrees of freedom. (Strictly speaking, 
we cannot speak of convergence rates since the algorithm will eventually fully resolve the 
problem.) 

Throughout this section we fix F — 1, N — 25000, and let be the Morse potential 

<f>(r) = exp(— 2a(r — 1)) — 2exp(— a(r — 1)), 

with the parameter a = 5. We defined the external force / to be 

-0.4^£ for -1/2 < x < 0, 

0.4^, for < x < 1/2. 

Note that f(x) behaves essentially like which is a typical decay rate for elastic fields 

generated by long-ranged defects in 2D/3D. (By contrast local perturbations decay exponen- 
tially in our ID model.) Thus, introducing this force allows us to study the performance of 
our adaptive algorithm in a setting that includes some of the features of 2D/3D problems. 
The constant 0.4 is fairly arbitrary. It was chosen sufficiently large to achieve a significant 
deformation, but sufficiently small to ensure that there exists an elastic stable equilibrium 
configuration. 

We will analyze the relative errors of the deformation gradient and of the energy defined, 
respectively, by 

IKc-yalU 2 (fi) , \E a (y a ) - E qc (y qc )\ 

and it-1/ \ r~i / f „m • (6-1) 



/(*) 



W*-F\\v<n) \ E M ~ E,(Fx 

In all computations, we compare the QC solutions against the (computable) exact solutions. 

6.1. A priori mesh refinement. We will compare the adaptive algorithm introduced in the 
next sub-section against a quasi-optimal a priori mesh refinement scheme, which exploits the 
known qualitative behaviour of the solution. For simplicity we keep the following discussions 
informal. 

We expect that, away from the defect, the exact solution will essentially behave like \u"\ < 
|/|. We choose to atomistic region the be of the form (—Me, Me). Closely following the 
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analysis of [2TJ Sec. 7.1] to optimize the mesh T h based on these assumptions, we obtain 
that the (quasi-)optimal mesh size in the continuum region is given, approximately, by 

f(Me) \x\ 



h*(x) = 2e 



f(x) Me ' 

The following algorithm generates a mesh with this quasi-optimal mesh size. We only state 
the construction for the right-hand half of the domain. The mesh will be symmetric about 
the centre x = 0. The factor 2e ensures that the restriction (T4) is observed. 

Algorithm 1 (A priori mesh refinement). 

(1) Add the vertices 0, s, . . . , eM to the mesh. 

(2) Let x\ be the right-most vertex already in the mesh. If x\ + h*(xfy > N, stop. 

(3) Otherwise, add the vertex x\ + h*{x\) to the mesh. Continue at (2). 

□ 

6.2. Adaptive algorithm. In our adaptive computations, we begin with a mesh that re- 
solves the "defect" (i.e., has atomistic resolution near x = 0) but is coarse elsewhere. We 
then employ the algorithm stated below to locally refine the mesh where required and thus 
improve the quality of the solution. 

Before we state the algorithm, we first define the error indicators according to which we 
drive the mesh refinement. Node-based error error constributions are split between neigh- 
bouring elements. Error contributions from the atomistic region are associated with the 
neighbouring continuum elements. 



The element error indicators for the gradient-driven algorithm are given by (cf. (3.7) and 



(3.13)) 



H-i + H + (vl) 2 + (vLi) 2 , if 4-i = R>, 

1 3^2 , 3„2 , /^f\2 



2 



Vk-i + hi + (vl) 2 , otherwise. 



Note that we have ignored the quadrature contributions. We have carefully justified this 
omission in Section 13.41 



The element error indicators for the energy-driven algorithm are given by (cf. Theorem 11 ) 



Mfe-i + \Pk + mLi + \A + fJ-l-i + t*l> if x k = R - 
Pk-=ASt?(Pk) 2 +{ Ivk-i+Pk + lpU + pl + pi + pU, ifx h k =L, 

\li k ~i + \nk + \Pk-i + \p{ + /-4> otherwise 



Since the energy error is formally second order, we did include the quadrature contribution 
\i q in this case. 

In the following algorithm, let p k G {Pk,pf }• Our algorithm is based on established ideas 
from the adaptive finite element literature [SI US] ■ 

Algorithm 2 (A posteriori mesh refinement). 

(1) Add the nodes 0, ±e, ...,±3e to the mesh. Add the nodes {x^ = —l/2,Xj C = 1/2} 
to the mesh. 

(2) Compute: Compute the QC solution on the current mesh, compute the error indi- 
cators pk- For T{l C f2 a , set pk = 0. 
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Figure 1. Relative errors in the deformation gradient (6.1) plotted against 



the number of degrees of freedom for three types of mesh refinements. 



(3) Mark: Choose a minimal subset Ai C {1, . . . , K} of indices such that 

1 K 

keM k=l 

(4) Refine: Bisect all elements T£ with indices belonging to A4. 

If an element that needs to be refined is adjacent to the atomistic region, merge this 
element into the atomistic region and create a new atomistic to continuum interface. 

(5) If the resulting mesh reaches a prescibed maximal number of degrees of freedom, stop 
algorithm; otherwise, go to Step (2). 

□ 

6.3. Numerical Results. We summarize the results of the computions with meshes gener- 
ated by the quasi-optimal refinement, and the adaptive algorithm with both gradient- and 
energy-based error indicators. 

(1) In Figure [I] we display the gradient errors for the three types of mesh generation 
algorithms: quasi-optimal a priori refinement, gradient-driven a posteriori refinement 
and energy-driven a posteriori refinement. The differences between the results pro- 
duced by the three algorithms is negligable. We observe rates close to (#DoF) _1 for 
all three algorithms. The efficiency indicators (estimate divided by actual error) are 
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10 2 10 3 10 4 

Number of Degrees of Freedom 



Figure 2. Efficiency factors (gradient a posteriori estimate divided by actual 
error) plotted against the number of degrees of freedom for three types of mesh 
refinements. 



displayed in Figure [2} They indicate an initial overestimation of the error, but after 
a sufficient accuracy is reached they enter a moderate range. 
(2) In Figure [3] we show the energy errors for the three types of mesh generation algo- 
rithms. Once again the differences between the three algorithms is negligable and 
the convergence rate is, as expected, twice the rate as compared with the gradient 
errors. However, the efficiency factors plotted in Figure [4] suggest that the constant 
prefactors in the estimators lead to a more substantial overestimation of the energy 
error. 

We can conclude that, at least in this model problem, both a posteriori error indicators 
can be used to select meshes that are quasi-optimal for both the deformation gradient and 
for the energy. In order to obtain sharper estimates on actual errors (in particular the energy 
error), alternative approaches such as the goal-oriented approach [31 [21] might be preferrable. 

7. Conclusion 

We derived a posteriori error estimators for the deformation gradient and for the energy in 
a ID atomistic chain, computed by an atomistic-to-continuum coupling method. Based on 
these estimators we proposed two adaptive mesh refinement algorithms, which we compared 
to quasi-optimal a priori mesh refinement. Our numerical experiments indicate that the 
resulting meshes are again quasi-optimal. 



A POSTERIORI ERROR CONTROL FOR A QC APPROXIMATION 



21 



10 u 



10" 



i 



10" 



10" 



10" 



10" 



10" 



A Priori 

Gradient-adaptivity 
Energy-adaptivity 




1.9393 



10 10 
Number of Degrees of Freedom 



10* 



Figure 3. Relative Error of the total energy (6.1 ) plotted against the number 



of degrees of freedom for three types of mesh refinements. 
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Figure 4. Efficiency Factor of the energy error estimator (energy a posteriori 
error estimate divided by actual energy error). 
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While further work in the application relevant 2D/3D setting remains to be done, we con- 
clude that adaptive mesh refinement driven by residual-based a posteriori error estimates can 
potentially lead to highly efficient atomistic/continuum multiscale computations of atomistic 
material defects. 
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